This invention relates generally to accelerating Monte Carlo calculations, and more particularly to determining radiation dose distributions in a patient for radiotherapy treatment planning by accelerating Monte Carlo dose distribution calculations through denoising of raw Monte Carlo dose distributions.
Medical equipment for radiation therapy treats tumorous (malignant) tissue with high energy radiation. The amount of radiation and its placement must be accurately controlled to ensure that the tumor receives sufficient radiation to be destroyed and that damage to the surrounding and adjacent non-tumorous tissue is minimized.
Internal source radiation therapy (referred to as brachytherapy) places capsules of radioactive material inside the patient in proximity to the tumorous tissue. Dose and placement are accurately controlled by the physical positioning of the isotope. However, internal source radiation therapy has disadvantages similar to those present with surgically invasive procedures, which include patient discomfort and risk of infection.
External source radiation therapy uses high energy radiation, collimated to direct a beam into the patient to the tumor site. Although the size and strength of the radiation beam from the external source may be accurately controlled outside of the patient, the dose received by a given volume within the patient may vary because of radiation scattering and absorption by intervening tissue. For these reasons, a determination of the proper dose and placement of the dose requires an estimation of the effects of treated tissue and the tissue surrounding the treated area in scattering or attenuation of the radiation beam.
The Monte Carlo method is the most accurate method for predicting dose distributions. The path of the particles (electrons, photons, protons, or neutrons) through the patient are simulated by Monte Carlo software. Particle interactions (scattering, attenuation, and energy deposition) are simulated one at a time. The Monte Carlo method traces paths of several million particles through a patient model, the patient model accurately reflecting the three dimensional variations of electron density within the volume of the patient under study. For large numbers of source radiation particles (typically above 107), the Monte Carlo method produces an accurate representation of the dose distribution. For these reasons, the Monte Carlo method is clinically preferred for the calculation of radiation dose in electron beam radiotherapy.
Unfortunately, the Monte Carlo method is extremely time consuming, taking upwards of an hour to compute a single dose distribution using current computers. One crucial and unique feature to Monte Carlo results is that the Monte Carlo method has no well-defined preset xe2x80x9cfinishxe2x80x9d time. Instead, dose distributions can be produced at any time, even after simulating only one source particle. However, the resulting dose distributions are xe2x80x9cnoisyxe2x80x9d (distorted), in inverse proportion to the number of source particles simulated. The Monte Carlo calculations are terminated when the noise level falls below a level deemed acceptable by the user.
Often, radiotherapy treatment planners wish to compare many dose distributions before selecting a final distribution for treatment. Hence, the need exists for a method of dose modeling which is as accurate as the Monte Carlo method but which has greater computational efficiency than the Monte Carlo method alone. This need especially exists as Monte Carlo programs are being developed for most commercial clinical systems.
More generally, the problem of noise is inherent in all Monte Carlo calculations that produce values on a grid. The problem of noise exists whether the grid elements are regularly spaced or irregularly spaced. For these reasons, there also exists a need to reduce noise and provide an accurate and computationally efficient estimate for Monte Carlo calculations in general. The present invention also satisfies this need, providing a method hereafter referred to as Monte Carlo denoising.
The present invention provides an estimate of an actual radiation dose distribution by reducing noise (denoising) from raw results of Monte Carlo generated dose distributions. The resulting denoised dose distribution more quickly and efficiently converges to the required level of accuracy for radiotherapy treatment planning than Monte Carlo calculations alone, and at a reduced computational cost. Optimal denoising reduces Monte Carlo run times by a factor of at least 3 to 6.
The present invention provides a method of accounting for details of the radiation source geometry and materials, and local changes in density of the patient at different points in the irradiated volume. This is done by denoising Monte Carlo dose distributions without significantly increasing the running time or introducing distortions to the underlying dose distributions.
The insight underlying the present invention is that a true (noiseless) radiation dose distribution is smoother (more spatially coherent) than a raw Monte Carlo result. In the present invention, it is the dose distribution itself that is accurately smoothed to reduce noise while not distorting the true underlying (signal) dose distribution.
It is our novel insight into the details of radiation transport physics which indicates that denoising is feasible and desirable (Deasy, June 2000, Denoising of electron beam Monte Carlo dose distributions using digital filtering techniques, Physics in Medicine and Biology, 45:1765-1779), which publication is incorporated in its entirety herein by reference. In simple terms, diffusive radiation transport implies that the actual dose distribution is smooth, whereas the Monte Carlo result without denoising contains statistical fluctuations which are more rough than the expected underlying dose distribution. Therefore, appropriate denoising techniques can be used to reduce the higher frequency noise of Monte Carlo results while not distorting the true underlying (signal) dose distribution.
In one aspect of the present invention, a method and an apparatus for accomplishing the method is provided where a raw Monte Carlo data distribution D0 (D0 everywhere refers to the raw Monte Carlo result) is first obtained and then denoised to produce a denoised distribution D1 (D1 everywhere refers to the denoised Monte Carlo result). This aspect of the present invention could be used to simply accelerate the computation of any Monte Carlo result or the present invention could be more specifically used to determine a dose distribution in a patient for radiation therapy treatment planning.
In another aspect of the present invention, denoising of the raw Monte Carlo dose distribution D0 could employ digital filtering, wavelet denoising, kernel smoothing or non-parametric regression smoothing. The digital filtering techniques could employ Binomial/Gaussian filters or local-least squares filters.
In another aspect of the present invention, denoising the raw Monte Carlo dose distribution D0 could first involve transforming the Monte Carlo dose distribution D0 by computing the square root of each data element of the Monte Carlo dose distribution D0. Digital filtering techniques could then be applied to the transformed data elements to obtain a filtered result. Then, the elements of the filtered result are squared to obtain a final best estimate of the dose distribution. The digital filtering techniques could employ Binomial/Gaussian filters or local-least squares filters. If employed, the Binomial/Gaussian filters could use a higher degree polynomial for less aggressive denoising. If employed, the local-least squares filters could use a smaller neighborhood for less aggressive denoising.
In a further aspect of the present invention, denoising the raw Monte Carlo dose distribution D0 employs either a spatially adaptive iterative filtering (SAIF) algorithm, a wavelet shrinkage threshold denoising algorithm, a spatially adaptive wavelet denoising (SAWD) algorithm, or a batch-averaged wavelet denoising (BAWD) algorithm.
In a still further aspect of the present invention, a method for determining a radiation dose pattern directed at a volume of interest in a patient includes providing to a computer a matrix T (describing an electron density for elements of the volume of interest in the patient), a subroutine (to produce a list of source particles describing a radiation source) and a Monte Carlo program (to simulate effects of the source particles). A dose accumulation matrix A is initialized to zeros. The computer then runs the Monte Carlo program to produce a dose deposition pattern for an incident source particle. The dose accumulation matrix A is then updated by adding the dose deposition pattern produced for the incident source particle. Running the Monte Carlo program and updating the dose accumulation matrix A is repeated until a predetermined stopping criteria is reached. The Monte Carlo dose distribution D0 is then set equal to the updated dose accumulation matrix A. Finally, a denoised estimate is developed for each element of the Monte Carlo data distribution D0 to produce a denoised dose distribution D1 representative of the radiation dose pattern. This aspect of the present invention could be varied by selecting an incident source particle to reduce uncertainty in a pre-determined location of a dose deposition pattern before running the Monte Carlo program. Selecting the incident source particle then becomes part of the repeated process (along with running the Monte Carlo program and updating the dose accumulation matrix A) until the pre-determined stopping criteria is reached.
In a yet further aspect of the present invention, developing the denoised estimate includes selecting an element from the Monte Carlo data distribution D0, creating a submatrix using values from the Monte Carlo data distribution D0 within a predetermined range of the selected element and multiplying the submatrix by a predetermined matrix designed to effect a least-squares fit then extracting the fit to the center element of the submatrix as the least-squares fit, and repeating these acts for each element of the Monte Carlo data distribution D0 to complete the development of the denoised estimate.
In another aspect of the present invention, developing the denoised estimate uses a variance stabilizing transformation (defined as a transformation which results in the statistical having more similar characteristics in every part of the data distribution). The variance stabilizing transformation could include computing the square root of each element of the Monte Carlo data distribution D0 to create a corresponding matrix D2. The corresponding matrix (D2) is then denoised to produce a matrix M. Every element of the matrix M is then squared to produce the denoised estimate.
In still another aspect of the present invention, developing the denoised estimate uses a statistically-based smoothing method. The statistically-based smoothing method could be the Loess method.
In yet another aspect of the present invention, developing the denoised estimate uses a curve fitting method employing Hermite polynomials for the least-squares fit. Developing the denoised estimate could also use a basis function selected from the group consisting of spline fits, spline fits with knots and weighted local least squares methods. Or, developing the denoised estimate could use a Fourier-based smoothing method. The Fourier-based smoothing method could employ a low-pass filter or a Wiener filter.
In a further aspect of the present invention, the Monte Carlo data distribution D0 is divided into parallel planes and developing the denoised estimate occurs independently for each plane.
In a still further aspect of the present invention, the Monte Carlo data distribution D0 is regridded onto a non-cartesian grid and developing the denoised estimate occurs on the non-cartesian grid with the resulting denoised dose distribution D1 being transformed back to a cartesian grid. The non-cartesian grid could be a diverging fan beam grid.
In a yet further aspect of the present invention, developing the denoised estimate uses a variable computational technique that varies to address selected portions of an image.
In another aspect of the present invention, developing the denoised estimate uses a computational technique that reduces random noise and accurately estimates underlying signal.
It is therefore one object of the invention to provide a method of dose calculation that is more computationally efficient at reaching the same level of noise than using present Monte Carlo methods alone to simulate radiation dose distributions.
It is another object of the invention to provide a method of dose calculation which reaches a lower level of dosage pattern noise than Monte Carlo methods for the same computation time.
It is a further object of the invention to accelerate any type of Monte Carlo calculation (for any purpose) to accurately represent a true underlying signal by generating and denoising a raw Monte Carlo calculation, the present invention reducing computation time relative to using present Monte Carlo methods alone.